The present invention relates to a semiconductor process device simulation method of simulating the manufacturing process and electrical characteristics of a semiconductor device using a computer, and a storage medium storing a simulation program and, more particularly, to a simulation method of efficiently and numerically solving, on a computer, simultaneous multi-dimensional linear equations based on physical laws used for a simulation.
As semiconductor integrated circuits such as a VLSI have become smaller, manufacturing has become increasingly complicated. As a result low-temperature processing, precise annealing processing, and accurate process design are of increasing importance.
On such a background, in recent years, device simulations are generally used for development and research of, e.g., the manufacture of semiconductor devices. A factor for this is a great progress in computers which actually execute the device simulation. Additionally, more convenient interfaces for the computers are becoming available. The progress in the analysis technique allows simulations for objects in a wider range and even analysis of breakdown phenomena or current concentration in a semiconductor device.
In computer simulations of thermal impurity diffusion which is one of the manufacturing processes of a semiconductor device described in reference 1: Ryo Dan, xe2x80x9cProcess Device Simulation Technologyxe2x80x9d, pp. 26-28, first a region to be analyzed is divided into meshes, and a diffusion equation is discretely defined for each mesh point. The diffusion equations are converted into linear equations and then into simultaneous linear equations by Newton""s method or the like so as to solve the diffusion equations.
The electrical characteristics of a semiconductor device prepared by manufacturing processes including thermal impurity diffusion, ion implantation, and thermal oxidation can be simulated by the method described in reference 2 (Dan Ryo, xe2x80x9cProcess Device Simulation Technologyxe2x80x9d, pp. 105-134), according to which a region to be analyzed is divided into meshes, and a Poisson equation and current continuous equation are discretely defined for each mesh point. These equations are converted into linear equations and then into simultaneous linear equations by Newton""s method or the like, and the equations are solved on a computer, thereby simulating the electrical characteristics of the semiconductor device.
In these semiconductor process device simulation methods, it is important in practice to solve large scale simultaneous linear equations having as many dimensions as the number of mesh points using a computer with as small a memory and as high a speed as possible.
The function of solving these simultaneous linear equations is called a xe2x80x9cmatrix solverxe2x80x9d for which various methods have been proposed. For example, for an equation having four dimensions factorization may be used.
And, for an equation having a large number of dimensions, an iterative method is used because of limitations on the memory capacity and calculation time needed. When the coefficient matrix of simultaneous linear equations to be solved is symmetrical, the ICCG (Incomplete Choleske and Conjugate Gradient) method is used. This is a CG (Conjugate-Gradient) method with preprocessing.
When the coefficient matrix is asymmetrical, a series of techniques called the Krylov subspace method are generally used because of advantages in convergence.
With both the ICCG and Krylov subspace methods, forming an algorithm for multiplying the original coefficient matrix with an approximate inverse matrix of the coefficient matrix of the equation to be solved permits a reduction in the number of conditions of the coefficient matrix, whereby the approximate inverse matrix of the coefficient matrix of the equation may be solved by a reduced number of iterations. This operation is called xe2x80x9cpreconditioningxe2x80x9d, as described above. The approximate inverse matrix of the coefficient matrix is called a xe2x80x9cpreconditioning matrixxe2x80x9d.
As the inverse matrix of the preconditioning matrix, a matrix obtained by xe2x80x9cincomplete LU-factorizationxe2x80x9d, i.e., LU-factorization, with limited fill-in generation positions is often used. Fill-in is a process of LU-factorization represented by equations (1) and (2) below and means that at a position (i, j) where an element value Aij is 0 in the original coefficient matrix, a new element Lij or Uij other than 0 is generated. When the generated fill-in is not rejected, and taken into consideration in the subsequent process of                               L                      i            ⁢                          xe2x80x83                        ⁢            j                          =                              A                          i              ⁢                              xe2x80x83                            ⁢              j                                -                                    ∑                              k                =                1                                            j                -                1                                      ⁢                          xe2x80x83                        ⁢                                          L                                  i                  ⁢                                      xe2x80x83                                    ⁢                  k                                            ⁢                              U                                  k                  ⁢                                      xe2x80x83                                    ⁢                  j                                            ⁢                              xe2x80x83                            ⁢                              (                                  k                   less than                   j                   less than                   i                                )                                                                        (        1        )                                          U                      i            ⁢                          xe2x80x83                        ⁢            j                          =                                            L                              i                ⁢                                  xe2x80x83                                ⁢                i                                            -                1                                      ⁡                          (                                                A                                      i                    ⁢                                          xe2x80x83                                        ⁢                    j                                                  -                                                      ∑                                          k                      =                      1                                                              i                      -                      1                                                        ⁢                                      xe2x80x83                                    ⁢                                                            L                                              i                        ⁢                                                  xe2x80x83                                                ⁢                        k                                                              ⁢                                          U                                              k                        ⁢                                                  xe2x80x83                                                ⁢                        j                                                                                                        )                                ⁢                      xe2x80x83                    ⁢                      (                          k               less than               i               less than               j                        )                                              (        2        )            
LU-factorization, conventional complete LU-factorization is performed.
In preprocessing based on incomplete LU-factorization, it is required to perform optimum processing which satisfies the following two contradictory conditions. First, fill-in generation should be minimized to suppress the amount of memory used and calculation time needed for performing the calculation. Second, a result as close as possible to complete LU-factorization should be obtained to increase convergence of the iterative method.
A conventional method of efficiently performing such processing is described in reference 3 (Shin Nakamura and Akio Nakagawa, High-speed Iterative Method for Two-dimensional Bipolar MOS Composite Device Simulator TONADDEIIxe2x80x9d, IEICE Technical Report, pp. 64-65) or reference 4 (Megumi Kawakami, Masahiro Sugaya, and Shiroo Kamohara, xe2x80x9cA New High-speed Non-equilibrium Point Defect Model for Annealing Simulationxe2x80x9d, SISPAD ""96, p. 94, FIG. 1a) is used.
In the method described in reference 3 or 4, as shown in the flow chart of FIG. 8, first, the maximum value n of the number of equations to be solved which are defined for each point is obtained in step 801. In step 802, the rows and columns of the coefficient matrix are put into a group in units of square submatrices of nxc3x97n for each mesh point, thereby forming blocks of coefficient matrices. In step 803, while virtually regarding an nxc3x97n square submatrix as one matrix, incomplete LU-factorization processing based on calculation of the square submatrices is performed. FIG. 9 shows the structure of the coefficient matrix formed by the above operation.
With the operation in units of square submatrices, incomplete LU-factorization close to complete LU-factorization can be more effectively performed using the same memory capacity as would be in processing units of scalar elements. This is because the fill-in generation pattern is determined while giving prominence to a combination of equations at one mesh point.
Objects to which the above-described conventional semiconductor process device simulation technology is applied are disclosed in Japanese Patent Laid-Open No. 6-53155 or 8-213334.
However, the above-described conventional semiconductor process device simulation technology has the following problem. That is, when a equation is derived from a certain physical law holds for only part of an analysis region to be processed, or no equation is defined, excess memory capacity is required to process incomplete LU-factorization.
This is because in making the computer calculations, a square submatrix having a uniform size of nxc3x97n is used as the processing unit of incomplete LU-factorization. For example, in the above-described simulation, the submatrix size of nxc3x97n is maintained even at a mesh point where n equations are not defined. The computer performing the calculations must form an overall coefficient matrix while inserting xe2x80x9c1xe2x80x9d to the corresponding principal diagonal portions of the submatrix and set the coefficient matrix on the memory. For this reason, excess memory capacity is used.
When the computer executes its calculation, excess memory is accessed. This wasteful processing adds to processing time.
It is the main object of the present invention to provide high speed simulation without using an excess memory for calculation.
It is another object of the present invention to provide a semiconductor process device simulation method which allows high-speed incomplete LU-factorization without using excess memory capacity even when an equation derived from a certain physical law and used to simulate the manufacturing process or electrical characteristics of a semiconductor device using a computer holds at only one portion of a region to be analyzed, or no equation is defined, and a storage medium storing a simulation program therefor.
In order to achieve the above objects, according to one aspect of the present invention there is provided a semiconductor device manufacturing process simulation method to aid manufacturers in forecasting electrical characteristics of semiconductor devices by performing a plurality of matrix manipulations of terms representing physical properties of the semiconductor devices, the matrices to be manipulated representing multidimensional simultaneous linear equations that are to be solved by a matrix solver that uses an iterative method in which the matrices are preconditioned by incomplete LU-factorization, the method for generating fill-ins for the matrices comprising the steps of: a first step of dividing the surface of a semiconductor device to be processed into a plurality of rectangles forming a mesh of predetermined size; a second step of assigning a numerical value to each mesh point of the mesh; a third step of setting equations representing a relationship among the plurality of numerical values; a fourth step of generating a coefficient matrix constituted by a plurality of principal diagonal submatrices each of which is arranged at each one of principal diagonal positions corresponding to each mesh point and representing a self feedback function to the mesh point, the coefficient matrix having rows and columns in numbers corresponding to the number of mesh points, and a plurality of non-principal diagonal submatrices each of which is arranged on any one of rows and columns passing through the principal diagonal positions corresponding to the mesh point and representing an interaction between the mesh point corresponding to the principal diagonal located on the same row or on the same column of the coefficient matrix with the non-principal diagonal submatrix and an adjacent mesh point connected to the mesh point through a mesh branch; a fifth step of performing calculation for the submatrices while regarding each submatrix of the coefficient matrix as one element, thereby performing incomplete LU-factorization of the coefficient matrix, wherein each of the principal diagonal submatrices is a square having rows and columns equal in number to equations set for a mesh point corresponding to the principal diagonal submatrix, each of the non-principal diagonal submatrices arranged in a row direction in correspondence with each of the mesh points is a matrix having rows equal in number to equations set at a mesh point corresponding to the principal diagonal submatrix located in the row and columns equal in number to equations set at an adjacent mesh point connected to the mesh point through a mesh branch, and each of the non-principal diagonal submatrices arranged in a column direction in correspondence with each of the mesh points is a matrix having columns equal in number to equations set at a mesh point corresponding to the principal diagonal submatrix located in the column and rows equal in number to the equations set at an adjacent mesh point connected to the mesh point through a mesh branch; and a sixth step of producing a signal indicative of the result of said calculation.
In another aspect of the invention there is provided a computer readable memory storing a semiconductor device manufacturing process simulation program to aid manufacturers in forecasting electrical characteristics of semiconductor devices by performing a plurality of matrix manipulations of terms representing physical properties of the semiconductor devices, the matrices to be manipulated representing multidimensional simultaneous linear equations that are to be solved by a matrix solver that uses an iterative method in which the matrices are preconditioned by incomplete LU-factorization, the program including a routine for generating fill-ins for the matrices, the routine comprising the steps of: causing a computer to perform the following functions: a first function of dividing a surface of a semiconductor device to be processed into a mesh of predetermined size; a second function of assigning a numerical value to each mesh point of the mesh; a third function of setting equations representing a relationship among the plurality of numerical values; a fourth function of generating a coefficient matrix constituted by a plurality of principal diagonal submatrices each of which is arranged at each one of principal diagonal positions corresponding to each mesh point and representing a self feedback function at the mesh point, the coefficient matrix having rows and columns in numbers corresponding to the number of mesh points, and a plurality of non-principal diagonal submatrices each of which is arranged on any one of rows and columns and representing an interaction between the mesh point corresponding to the principal diagonal positions located on the same row or on the same column of the coefficient matrix with the non-principle diagonal submatrix and an adjacent mesh point connected to the mesh point through a mesh branch; and a fifth function of performing calculation for the submatrices while regarding each submatrix of the coefficient matrix as one element; thereby performing incomplete LU-factorization of the coefficient matrix, wherein each of the principal diagonal submatrices is a square matrix having rows and columns equal in number to equations set forth for a mesh point corresponding to the principal diagonal submatrix, each of the non-principal diagonal submatrices arranged in a row direction in correspondence with each of the mesh points is a matrix having rows equal in number to equations set at a mesh point corresponding to the principal diagonal submatrix and located in the row and columns equal in number to equations set at an adjacent mesh point connected to the mesh point through a mesh branch, and each of the non-principal diagonal submatrices arranged in a column direction in correspondence with each of the mesh points is a matrix having columns equal in number to equations set a mesh point corresponding tot he principal diagonal submatrix located in the column and rows equal in number to equations set at an adjacent mesh point connected to the mesh point through a mesh branch.
With this arrangement, the non-principal diagonal submatrix not always will be a square matrix. When the number of equations set at an adjacent mesh point connected to a mesh point of interest through a mesh branch is smaller than that set for the mesh point of interest, a non-square matrix is obtained, so the number of calculations is reduced as compared to a case in which a square matrix is simply assigned.